# PoP - Policing Socio-Geographic Boundaries and Inequality
# script for creating Figure 6
# Figure 6 - Heterogeneous Influence of Racial Boundary on Police Stops By Race of Civilian. 


# load packages
suppressPackageStartupMessages(
  
  {
    library(dplyr)
    library(tidyverse)
    library(ggplot2)
    library(haven)
    library(readxl)
    library(readr)
    library(areal)
    library(car)
    library(estimatr)
    library(magrittr)
    library(texreg)
    library(sandwich)
    library(jtools)
    library(ggh4x)
  }
)

## Austin

load("aus_stops_final.RData")

# log and +1 to variables (pop already logged)
aus_stops_final = aus_stops_final %>% 
  mutate(lmhhi = log(mhhi + 1),
         larrests = log(total_arrests + 1),
         lmisdemeanors = log(misdemeanor_arrests + 1),
         lfelonies = log(felony_arrests + 1),
         lnonviolent = log(nonviolent_arrests + 1),
         lviolent = log(violent_arrests + 1),
         lsociety = log(society_arrests + 1),
         lperson = log(person_arrests + 1),
         lproperty = log(property_arrests + 1),
         lcrime = log(crime_all + 1),
         lpropertycrime = log(crime_property + 1),
         lviolentcrime = log(crime_violent + 1))

# now for new stop vars
aus_stops_final = aus_stops_final %>%
  mutate(lall_stops_total = log(all_stops_total + 1),
         lall_stops_black = log(all_stops_black + 1),
         lall_stops_latino = log(all_stops_latino + 1),
         lall_stops_white = log(all_stops_white + 1),
         lall_stops_asian = log(all_stops_asian + 1),
         lall_stops_nonwhite = log(all_stops_nonwhite + 1))

# create scaled DVs (to account for different pop. sizes and density in blocks)
aus_stops_final$larrest_sd <- aus_stops_final$larrests - (mean(aus_stops_final$larrests)/sd(aus_stops_final$larrests))
aus_stops_final$lmisdemeanors_sd <- aus_stops_final$lmisdemeanors - (mean(aus_stops_final$lmisdemeanors)/sd(aus_stops_final$lmisdemeanors))
aus_stops_final$lfelonies_sd <- aus_stops_final$lfelonies - (mean(aus_stops_final$lfelonies)/sd(aus_stops_final$lfelonies))
aus_stops_final$lnonviolent_sd <- aus_stops_final$lnonviolent - (mean(aus_stops_final$lnonviolent)/sd(aus_stops_final$lnonviolent))
aus_stops_final$lviolent_sd <- aus_stops_final$lviolent - (mean(aus_stops_final$lviolent)/sd(aus_stops_final$lviolent))
aus_stops_final$lsociety_sd <- aus_stops_final$lsociety - (mean(aus_stops_final$lsociety)/sd(aus_stops_final$lsociety))
aus_stops_final$lperson_sd <- aus_stops_final$lperson - (mean(aus_stops_final$lperson)/sd(aus_stops_final$lperson))
aus_stops_final$lproperty_sd <- aus_stops_final$lproperty - (mean(aus_stops_final$lproperty)/sd(aus_stops_final$lproperty))

# create scaled DVs for stop vars
aus_stops_final$lall_stops_total_sd <- aus_stops_final$lall_stops_total - (mean(aus_stops_final$lall_stops_total)/sd(aus_stops_final$lall_stops_total))
aus_stops_final$lall_stops_black_sd <- aus_stops_final$lall_stops_black - (mean(aus_stops_final$lall_stops_black)/sd(aus_stops_final$lall_stops_black))
aus_stops_final$lall_stops_latino_sd <- aus_stops_final$lall_stops_latino - (mean(aus_stops_final$lall_stops_latino)/sd(aus_stops_final$lall_stops_latino))
aus_stops_final$lall_stops_white_sd <- aus_stops_final$lall_stops_white - (mean(aus_stops_final$lall_stops_white)/sd(aus_stops_final$lall_stops_white))
aus_stops_final$lall_stops_asian_sd <- aus_stops_final$lall_stops_asian - (mean(aus_stops_final$lall_stops_asian)/sd(aus_stops_final$lall_stops_asian))
aus_stops_final$lall_stops_nonwhite_sd <- aus_stops_final$lall_stops_nonwhite - (mean(aus_stops_final$lall_stops_nonwhite)/sd(aus_stops_final$lall_stops_nonwhite))


# create binned boundary measure
aus_stops_final <- aus_stops_final %>% mutate(boundary_quart = quantile(aus_stops_final$p_race_white_blv, prob=.75, na.rm=TRUE),
                                              boundary_quart_dummy = ifelse(p_race_white_blv > boundary_quart,1,0 ))

# Chicago

load("chi_stops_final.RData")

# now create binned racial blv measures
chi_stops_final <- chi_stops_final %>% mutate(boundary_quart = quantile(chi_stops_final$p_race_white_blv, prob=.75, na.rm=TRUE),
                                              boundary_quart_dummy = ifelse(p_race_white_blv > boundary_quart,1,0 ))


# log and +1 to variables (pop already logged)
chi_stops_final = chi_stops_final %>% 
  mutate(lmhhi = log(mhhi + 1),
         larrests = log(total_arrests + 1),
         lmisdemeanors = log(misdemeanor_arrests + 1),
         lfelonies = log(felony_arrests + 1),
         lnonviolent = log(nonviolent_arrests + 1),
         lviolent = log(violent_arrests + 1),
         lsociety = log(society_arrests + 1),
         lperson = log(person_arrests + 1),
         lproperty = log(property_arrests + 1),
         lcrime = log(crime_all + 1),
         lpropertycrime = log(property_crime + 1),
         lviolentcrime = log(violent_crime + 1))

# now for new stop vars
chi_stops_final = chi_stops_final %>%
  mutate(lall_stops_total = log(all_stops_total + 1),
         lall_stops_black = log(all_stops_black + 1),
         lall_stops_latino = log(all_stops_latino + 1),
         lall_stops_white = log(all_stops_white + 1),
         lall_stops_asian = log(all_stops_asian + 1),
         lall_stops_nonwhite = log(all_stops_nonwhite + 1))

# create scaled DVs (to account for different pop. sizes and density in blocks)
chi_stops_final$larrest_sd <- chi_stops_final$larrests - (mean(chi_stops_final$larrests)/sd(chi_stops_final$larrests))
chi_stops_final$lmisdemeanors_sd <- chi_stops_final$lmisdemeanors - (mean(chi_stops_final$lmisdemeanors)/sd(chi_stops_final$lmisdemeanors))
chi_stops_final$lfelonies_sd <- chi_stops_final$lfelonies - (mean(chi_stops_final$lfelonies)/sd(chi_stops_final$lfelonies))
chi_stops_final$lnonviolent_sd <- chi_stops_final$lnonviolent - (mean(chi_stops_final$lnonviolent)/sd(chi_stops_final$lnonviolent))
chi_stops_final$lviolent_sd <- chi_stops_final$lviolent - (mean(chi_stops_final$lviolent)/sd(chi_stops_final$lviolent))
chi_stops_final$lsociety_sd <- chi_stops_final$lsociety - (mean(chi_stops_final$lsociety)/sd(chi_stops_final$lsociety))
chi_stops_final$lperson_sd <- chi_stops_final$lperson - (mean(chi_stops_final$lperson)/sd(chi_stops_final$lperson))
chi_stops_final$lproperty_sd <- chi_stops_final$lproperty - (mean(chi_stops_final$lproperty)/sd(chi_stops_final$lproperty))

# create scaled DVs for stop vars
chi_stops_final$lall_stops_total_sd <- chi_stops_final$lall_stops_total - (mean(chi_stops_final$lall_stops_total)/sd(chi_stops_final$lall_stops_total))
chi_stops_final$lall_stops_black_sd <- chi_stops_final$lall_stops_black - (mean(chi_stops_final$lall_stops_black)/sd(chi_stops_final$lall_stops_black))
chi_stops_final$lall_stops_latino_sd <- chi_stops_final$lall_stops_latino - (mean(chi_stops_final$lall_stops_latino)/sd(chi_stops_final$lall_stops_latino))
chi_stops_final$lall_stops_white_sd <- chi_stops_final$lall_stops_white - (mean(chi_stops_final$lall_stops_white)/sd(chi_stops_final$lall_stops_white))
chi_stops_final$lall_stops_asian_sd <- chi_stops_final$lall_stops_asian - (mean(chi_stops_final$lall_stops_asian)/sd(chi_stops_final$lall_stops_asian))
chi_stops_final$lall_stops_nonwhite_sd <- chi_stops_final$lall_stops_nonwhite - (mean(chi_stops_final$lall_stops_nonwhite)/sd(chi_stops_final$lall_stops_nonwhite))

## Milwaukee
load("mil_stops_final.RData")

# log and +1 to variables (pop already logged)
mil_stops_final = mil_stops_final %>% 
  mutate(lmhhi = log(mhhi + 1),
         larrests = log(total_arrests + 1),
         lmisdemeanors = log(misdemeanor_arrests + 1),
         lfelonies = log(felony_arrests + 1),
         lnonviolent = log(nonviolent_arrests + 1),
         lviolent = log(violent_arrests + 1),
         lsociety = log(society_arrests + 1),
         lperson = log(person_arrests + 1),
         lproperty = log(property_arrests + 1),
         lcrime = log(crime_all + 1),
         lpropertycrime = log(crime_property + 1),
         lviolentcrime = log(crime_violent + 1))

# now for new stop vars
mil_stops_final = mil_stops_final %>%
  mutate(lall_stops_total = log(all_stops_total + 1),
         lall_stops_black = log(all_stops_black + 1),
         lall_stops_latino = log(all_stops_latino + 1),
         lall_stops_white = log(all_stops_white + 1),
         lall_stops_asian = log(all_stops_asian + 1),
         lall_stops_nonwhite = log(all_stops_nonwhite + 1),
         lped_stops_total = log(ped_stops_total + 1),
         lped_stops_black = log(ped_stops_black + 1),
         lped_stops_latino = log(ped_stops_latino + 1),
         lped_stops_white = log(ped_stops_white + 1),
         lped_stops_asian = log(ped_stops_asian + 1),
         lped_stops_nonwhite = log(ped_stops_nonwhite + 1))

# create scaled DVs (to account for different pop. sizes and density in blocks)
mil_stops_final$larrest_sd <- mil_stops_final$larrests - (mean(mil_stops_final$larrests)/sd(mil_stops_final$larrests))
mil_stops_final$lmisdemeanors_sd <- mil_stops_final$lmisdemeanors - (mean(mil_stops_final$lmisdemeanors)/sd(mil_stops_final$lmisdemeanors))
mil_stops_final$lfelonies_sd <- mil_stops_final$lfelonies - (mean(mil_stops_final$lfelonies)/sd(mil_stops_final$lfelonies))
mil_stops_final$lnonviolent_sd <- mil_stops_final$lnonviolent - (mean(mil_stops_final$lnonviolent)/sd(mil_stops_final$lnonviolent))
mil_stops_final$lviolent_sd <- mil_stops_final$lviolent - (mean(mil_stops_final$lviolent)/sd(mil_stops_final$lviolent))
mil_stops_final$lsociety_sd <- mil_stops_final$lsociety - (mean(mil_stops_final$lsociety)/sd(mil_stops_final$lsociety))
mil_stops_final$lperson_sd <- mil_stops_final$lperson - (mean(mil_stops_final$lperson)/sd(mil_stops_final$lperson))
mil_stops_final$lproperty_sd <- mil_stops_final$lproperty - (mean(mil_stops_final$lproperty)/sd(mil_stops_final$lproperty))

# create scaled DVs for stop vars
mil_stops_final$lall_stops_total_sd <- mil_stops_final$lall_stops_total - (mean(mil_stops_final$lall_stops_total)/sd(mil_stops_final$lall_stops_total))
mil_stops_final$lall_stops_black_sd <- mil_stops_final$lall_stops_black - (mean(mil_stops_final$lall_stops_black)/sd(mil_stops_final$lall_stops_black))
mil_stops_final$lall_stops_latino_sd <- mil_stops_final$lall_stops_latino - (mean(mil_stops_final$lall_stops_latino)/sd(mil_stops_final$lall_stops_latino))
mil_stops_final$lall_stops_white_sd <- mil_stops_final$lall_stops_white - (mean(mil_stops_final$lall_stops_white)/sd(mil_stops_final$lall_stops_white))
mil_stops_final$lall_stops_asian_sd <- mil_stops_final$lall_stops_asian - (mean(mil_stops_final$lall_stops_asian)/sd(mil_stops_final$lall_stops_asian))
mil_stops_final$lall_stops_nonwhite_sd <- mil_stops_final$lall_stops_nonwhite - (mean(mil_stops_final$lall_stops_nonwhite)/sd(mil_stops_final$lall_stops_nonwhite))


# now create binned racial blv measures
mil_stops_final <- mil_stops_final %>% mutate(boundary_quart = quantile(mil_stops_final$p_race_white_blv, prob=.75, na.rm=TRUE),
                                              boundary_quart_dummy = ifelse(p_race_white_blv > boundary_quart,1,0 ))


###### race_blvXp_white plots w race of person stopped #####
## Austin

form_list = 
  list(
    'lall_stops_nonwhite_sd ~ boundary_quart_dummy*p_race_white + ses_blv + lpop + age_15_35_male + hhi + lmhhi + 
    pown + p_poverty + p_emp_unemployed + pcol + lpropertycrime + lviolentcrime + com_density + phys_bound',
    'lall_stops_white_sd ~ boundary_quart_dummy*p_race_white + ses_blv + lpop + age_15_35_male + hhi + lmhhi + 
    pown + p_poverty + p_emp_unemployed + pcol + lpropertycrime + lviolentcrime + com_density + phys_bound'
  )

hetdf_list = list(form_list)
hetdf_lab = c("Non-white Persons","White Persons")

out_mod_list = as.list(rep(NA, length(form_list)))
out_df_list = as.list(rep(NA, length(form_list)))

for (i in 1:length(form_list)) {
  
  print(paste0("Iteration ", i))
  
  out_mod_list[[i]] = lm_robust(as.formula(form_list[[i]]), aus_stops_final)
  
  fakedata = 
    aus_stops_final[, names(out_mod_list[[i]]$coefficients)[2:(length(names(out_mod_list[[i]]$coefficients))-1)]] %>% 
    filter(lpop > 0) %>% 
    as.data.frame %>% 
    mutate(geometry = NULL) %>% 
    apply(X = ., MARGIN = 2, FUN = function(x) mean(x, na.rm = TRUE)) %>% 
    t %>% 
    as.data.frame %>% 
    slice(rep(1:n(), each = 22)) %>% 
    mutate(p_race_white = rep(seq(from = 0, to = 1, by = .1), 2),
           boundary_quart_dummy = c(rep(0, 11), rep(1, 11)))
  
  out_df_list[[i]] = data.frame(
    
    est = predict(object = out_mod_list[[i]], newdata = fakedata, se.fit = TRUE)$fit,
    se = predict(object = out_mod_list[[i]], newdata = fakedata, se.fit = TRUE)$se.fit,
    p_race_white = rep(seq(from = 0, to = 1, by = .1), 2),
    boundary_quart_dummy = c(rep(0, 11), rep(1, 11))
    
  ) %>% 
    mutate(dataset = hetdf_lab[i])
  
}

df_to_plot = out_df_list %>%
  do.call(rbind.data.frame, .) %>%
  mutate(boundary_quart_dummy = factor(boundary_quart_dummy, labels = c("No Boundary", "Boundary")))

df_to_plot = df_to_plot %>%
  mutate(city = "Austin")

# create plot
aus_stops_fow_plot = df_to_plot %>%
  ggplot() +
  geom_point(aes(x = p_race_white, y = est, col = boundary_quart_dummy),
             position = position_dodge(.075),
             size = 2) +
  geom_errorbar(aes(x = p_race_white,
                    ymin = est - 1.96 * se,
                    ymax = est + 1.96 * se,
                    col = boundary_quart_dummy),
                width = 0,
                size = .4,
                position = position_dodge(.075)) +
  geom_errorbar(aes(x = p_race_white,
                    ymin = est - 1.645 * se,
                    ymax = est + 1.645 * se,
                    col = boundary_quart_dummy),
                width = 0,
                size = .6,
                position = position_dodge(.075)) +
  scale_color_grey(start = .6, end = 0) +
  labs(x = "% White",
       y = "Predicted Value",
       col = "Race Boundary") +
  facet_grid2(city ~ dataset, scales = "free_x", independent = "x") +
  theme_bw(base_size=12,base_family="Times") +
  theme(panel.border=element_rect(fill=NA, colour=NA), 
        legend.title = element_blank(),
        legend.position="none",
        panel.grid.minor=element_line(colour=NA)) +
  theme(strip.background = element_rect(fill="white")) + 
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        axis.title.x=element_blank())


## Chicago
form_list = 
  list(
    'lall_stops_nonwhite_sd ~ boundary_quart_dummy*p_race_white + ses_blv + lpop + age_15_35_male + hhi + lmhhi + 
    pown + p_poverty + p_emp_unemployed + pcol + lpropertycrime + lviolentcrime + com_density + phys_bound',
    'lall_stops_white_sd ~ boundary_quart_dummy*p_race_white + ses_blv + lpop + age_15_35_male + hhi + lmhhi + 
    pown + p_poverty + p_emp_unemployed + pcol + lpropertycrime + lviolentcrime + com_density + phys_bound'
  )

hetdf_list = list(form_list)
hetdf_lab = c("Non-white Persons","White Persons")

out_mod_list = as.list(rep(NA, length(form_list)))
out_df_list = as.list(rep(NA, length(form_list)))

for (i in 1:length(form_list)) {
  
  print(paste0("Iteration ", i))
  
  out_mod_list[[i]] = lm_robust(as.formula(form_list[[i]]), chi_stops_final)
  
  fakedata = 
    chi_stops_final[, names(out_mod_list[[i]]$coefficients)[2:(length(names(out_mod_list[[i]]$coefficients))-1)]] %>% 
    filter(lpop > 0) %>% 
    as.data.frame %>% 
    mutate(geometry = NULL) %>% 
    apply(X = ., MARGIN = 2, FUN = function(x) mean(x, na.rm = TRUE)) %>% 
    t %>% 
    as.data.frame %>% 
    slice(rep(1:n(), each = 22)) %>% 
    mutate(p_race_white = rep(seq(from = 0, to = 1, by = .1), 2),
           boundary_quart_dummy = c(rep(0, 11), rep(1, 11)))
  
  out_df_list[[i]] = data.frame(
    
    est = predict(object = out_mod_list[[i]], newdata = fakedata, se.fit = TRUE)$fit,
    se = predict(object = out_mod_list[[i]], newdata = fakedata, se.fit = TRUE)$se.fit,
    p_race_white = rep(seq(from = 0, to = 1, by = .1), 2),
    boundary_quart_dummy = c(rep(0, 11), rep(1, 11))
    
  ) %>% 
    mutate(dataset = hetdf_lab[i])
  
}

df_to_plot = out_df_list %>%
  do.call(rbind.data.frame, .) %>%
  mutate(boundary_quart_dummy = factor(boundary_quart_dummy, labels = c("No Boundary", "Boundary")))

df_to_plot = df_to_plot %>%
  mutate(city = "Chicago")

# create plot
chi_stops_fow_plot = df_to_plot %>%
  ggplot() +
  geom_point(aes(x = p_race_white, y = est, col = boundary_quart_dummy),
             position = position_dodge(.075),
             size = 2) +
  geom_errorbar(aes(x = p_race_white,
                    ymin = est - 1.96 * se,
                    ymax = est + 1.96 * se,
                    col = boundary_quart_dummy),
                width = 0,
                size = .4,
                position = position_dodge(.075)) +
  geom_errorbar(aes(x = p_race_white,
                    ymin = est - 1.645 * se,
                    ymax = est + 1.645 * se,
                    col = boundary_quart_dummy),
                width = 0,
                size = .6,
                position = position_dodge(.075)) +
  scale_color_grey(start = .6, end = 0) +
  labs(x = "% White",
       y = "Predicted Value",
       col = "Race Boundary") +
  facet_grid2(city ~ dataset, scales = "free_x", independent = "x") +
  theme_bw(base_size=12,base_family="Times") +
  theme(panel.border=element_rect(fill=NA, colour=NA), 
        legend.title = element_blank(),
        legend.position="none", # remove legend
        panel.grid.minor=element_line(colour=NA)) +
  theme(strip.background = element_rect(fill="white")) + # remove grey from background
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), # remove gridlines from background
        axis.title.x=element_blank(), strip.text.x = element_blank()) # remove x axis title and x labels for plots




## Milwaukee
form_list = 
  list(
    'lall_stops_nonwhite_sd ~ boundary_quart_dummy*p_race_white + ses_blv + lpop + age_15_35_male + hhi + lmhhi + 
    pown + p_poverty + p_emp_unemployed + pcol + lpropertycrime + lviolentcrime + com_density + phys_bound',
    'lall_stops_white_sd ~ boundary_quart_dummy*p_race_white + ses_blv + lpop + age_15_35_male + hhi + lmhhi + 
    pown + p_poverty + p_emp_unemployed + pcol + lpropertycrime + lviolentcrime + com_density + phys_bound'
  )

hetdf_list = list(form_list)
hetdf_lab = c("Non-white Persons","White Persons")

out_mod_list = as.list(rep(NA, length(form_list)))
out_df_list = as.list(rep(NA, length(form_list)))

for (i in 1:length(form_list)) {
  
  print(paste0("Iteration ", i))
  
  out_mod_list[[i]] = lm_robust(as.formula(form_list[[i]]), mil_stops_final)
  
  fakedata = 
    mil_stops_final[, names(out_mod_list[[i]]$coefficients)[2:(length(names(out_mod_list[[i]]$coefficients))-1)]] %>% 
    filter(lpop > 0) %>% 
    as.data.frame %>% 
    mutate(geometry = NULL) %>% 
    apply(X = ., MARGIN = 2, FUN = function(x) mean(x, na.rm = TRUE)) %>% 
    t %>% 
    as.data.frame %>% 
    slice(rep(1:n(), each = 22)) %>% 
    mutate(p_race_white = rep(seq(from = 0, to = 1, by = .1), 2),
           boundary_quart_dummy = c(rep(0, 11), rep(1, 11)))
  
  out_df_list[[i]] = data.frame(
    
    est = predict(object = out_mod_list[[i]], newdata = fakedata, se.fit = TRUE)$fit,
    se = predict(object = out_mod_list[[i]], newdata = fakedata, se.fit = TRUE)$se.fit,
    p_race_white = rep(seq(from = 0, to = 1, by = .1), 2),
    boundary_quart_dummy = c(rep(0, 11), rep(1, 11))
    
  ) %>% 
    mutate(dataset = hetdf_lab[i])
  
}

df_to_plot = out_df_list %>%
  do.call(rbind.data.frame, .) %>%
  mutate(boundary_quart_dummy = factor(boundary_quart_dummy, labels = c("No Boundary", "Boundary")))

df_to_plot = df_to_plot %>%
  mutate(city = "Milwaukee")

# create plot
mil_stops_fow_plot = df_to_plot %>%
  ggplot() +
  geom_point(aes(x = p_race_white, y = est, col = boundary_quart_dummy),
             position = position_dodge(.075),
             size = 2) +
  geom_errorbar(aes(x = p_race_white,
                    ymin = est - 1.96 * se,
                    ymax = est + 1.96 * se,
                    col = boundary_quart_dummy),
                width = 0,
                size = .4,
                position = position_dodge(.075)) +
  geom_errorbar(aes(x = p_race_white,
                    ymin = est - 1.645 * se,
                    ymax = est + 1.645 * se,
                    col = boundary_quart_dummy),
                width = 0,
                size = .6,
                position = position_dodge(.075)) +
  scale_color_grey(start = .6, end = 0) +
  labs(x = "% White",
       y = "Predicted Value",
       col = "Race Boundary") +
  facet_grid2(city ~ dataset, scales = "free_x", independent = "x") +
  theme_bw(base_size=12,base_family="Times") +
  theme(panel.border=element_rect(fill=NA, colour=NA), 
        legend.title = element_blank(),
        legend.position="bottom",
        panel.grid.minor=element_line(colour=NA)) +
  theme(strip.background = element_rect(fill="white")) + 
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        strip.text.x = element_blank())


# put all the plots together
library(cowplot)
stops_fow_plots <- plot_grid(aus_stops_fow_plot, chi_stops_fow_plot, mil_stops_fow_plot,
                             align = "v", nrow = 3, rel_heights = c(1/3, 1/3, 1/2.5))

ggsave(plot = stops_fow_plots, filename = "figure6_plots.png", 
       width = 8, height = 12, bg = 'white')


